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three parts. In part I we discuss some basic ingredients from the spectral Fourier and Cheby- 
shev approximation theory. Part II contains a brief survey on hyperbolic and parabolic time- 
dependent problems which are dealt with both the energymethod and the related Fourier 
analysis. In part III we combine the ideas presented in the first two parts, in our study of 
accuracy stability and convergence of the spectral Fourier approximation to time-dependent 
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1. SPECTRAL APPROXIMATION THEORY 


1.1. The Periodic Problem - Fourier Approximation 

Consider the first order Sturm-Liouville (SL) problem 
(1.1.1) -^<f> = A$(x), 0 < x < 2ir, 

augmented by periodic boundary conditions 


( 1 . 1 . 2 ) 


m = <t>{ 2 tt). 


It has an infinite sequence of eigenvalues, A* = ik } with the corresponding eigenfunctions 
<f>k(x) = e xkx . Thus, (A*. = ik , = e lkx ) are the eigenpairs of the differentiation operator 

D = ■£; in L 2 [0,27r), and they form a complete system in this space. This system is complete 
in the following sense. Let L 2 [0,27 t) be induced with the usual Euclidean inner product 

(1.1.3) (u^x),^®)) = / wi(x)w 2 (x)dx. 

Jo 

Note that <f>k{%) are orthogonal with respect to this inner product, for 

t 11 - 4 ) = = ]t k k. 

Let w{x)eL 2 [ 0,27t) be associated with its spectral representation in this system, i.e., the 
Fourier expansion 

(1.1.5) w(x ) ~ E w(k)<f> k (x), w(k) = 

k=-oo llwlr 

or equivalently, 

°° 1 r 2 * 

(1.1.6) w(x ) ~ E w(k)e tkx , w(k ) = — / tu(£)e -1 **. 

*=-oo 2v J t=° 


The truncated Fourier expansion 

(1.1.7) E w(k)e ikx , 

k=-N 
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denotes the spectral-Fourier projection of w(x) into 7r^-the space of trigonometric polyno- 
mials of degree < N : 


N 

s N w = *(o) + J2l™( k ) eikx + ™(- k ) e ~ ikx } 

Jb=l 

N 

(1.1-8) = u’(0) + + w(—k)] cos kx + — w( 

k=i 

N 

= ^ f (ik cos kx + bk sin kx\ 

k=o 

here a*, and b * are the usual Fourier coefficients given by 

1 /*2ir 

a k = u/(fc) 4 u;( — /c) = — w (0 cos k£d(, 

7T Vo 

(1.1.9) 

6* = i[u>(Ai) — w(— A:)] = — f w(£)smk£d£. 

i r Jo 

Since w — S^w is orthogonal to the 7Tjv-space: 

(1.1.10) (u; — e tkx ) = 27ru;(fc) — 2irw(k) = 0, |fc| < iV, 
it follows that for any pn^n we have (see Figure 1) 

(1.1.11) ||u> - pjv|| 2 = ||u> - 5 atu>|| 2 + 115^ - p w || 2 . 

Hence, solves the least-squares problem 


w 



(1.1.12) \\w - 5jvu;|| - Min PArCTjv ||iw - p w || 


A:)] sin kx 
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i.e., S N w is the best least-squares approximation to w. Moreover, (1.1.11) with p N - 0 yields 

(1.1.13) \\Snw\\ 2 = IMI 2 - lb - SjvHI 2 < HI 2 

and by letting N — > oo we arrive at Bessel’s inequality 

(U.14) 2 * e i*(*)i ! = E w*)l’IIM’ < IMI ! . 

Jb=-oo k=-oo 

Remark : An immediate consequence of (1.1.14) is the Riemann-Lebesgue lemma, asserting 

that 2ir 

w(k) = — / w(()e~ tk *d ( — > 0, for any weL 2 [0, 27r). 

2lT Jo k-+ oo 

The system {<f> k = e* fcl } is complete in the sense that for any w(x)eL 2 [0, 2tt) we have 
Parseval equality: 

(1.1.15) 2* E |*(*)| J = E MJOI’MI’ = IMIt 

Jk=-oo Jc=~oo 

which in view of (1.1.13), is the same as 

(1.1.16) lim || S N w = Y^w{k)e ikx - in(x)|| = 0. 

N-.00 _ N 

The last equality establishes the L 2 convergence of the spectral-Fourier projection, Snw(x), 
to u;(x), whose difference can be upper bounded by the following 

Error Estimate : 

l|io — 5Vtn|| 2 = ||H| 3 — II *s’jv' u; 1I 2 = Y* b(*)l 2 M 2 = ^ 'Jh b( fc )l • 

|fc|>Ar \k\>N 

We observe that the RHS tends to zero as a tail of a converging sequence, i.e., 

(1.1.17) \w(x) — Y w(k)e ikx \ 2 dx = 2tt Y H*)! 2 ^ 0 - 

Jo k=-N \k\>N N ^°° 

The last equality tells us that the convergence rate depends on how fast the Fourier coeffi- 
cients, w(k), decay to zero, and we shall quantify this in a more precise way below. 

Remark (1.1.17) yields uniform a.e., convergence for subsequences ; in fact one can show 

(1.1.18) a.e. ton |to(x) - S Wj) u;(x)| = 0, inf^-^l. 
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In fact, w(x ) — a.e. Iim7v_,oo Sprw(x) for all iweX 2 [0, 2tt], but a.e. convergence may fail if u;(-) 
is only L 1 [0, 27r]-integrable. 

Yet, if we agree to assume sufficient smoothness, we find the convergence of spectral- 
Fourier projection to be very rapid, both in the L 2 and the pointwise sense. To this we 
proceed as follows. 


Define the Sobolev space H*[0, 2ir) consisting of 27r-periodic functions which their first 
a-derivatives are T 2 -integrable; set the inner product 


(1.1.19) 


* y*2x 

(wi,w 2 )h‘ ~ z2 J D p Wi(x)D p w 2 (x)da 

p—Q J ® 


The essential ingredient here is that the system {e ,fcx } - which was already shown to be 
complete in L 2 [0,27 t) = ff°[0,27r), is also a complete system in H'[0,2n) for any a > 0. For 
orthogonality we have 


( 1 . 1 . 20 ) 




0 j ^ k, 
2tt Y^k 2p j = k. 

P= 0 


The Fourier expansion now reads 

(1.1.21) w(x) ~ £ w.(k)e ikx 


k——c 


where the Fourier coefficients, are given by 

(u/(a:), e***)^. 


( 1 . 1 . 22 ) 


w t (k) = 


(e ik *,e ikx ) H . ‘ 
We integrate by parts and use periodicity to obtain 

* f2w 


(tw(x), e ,fcl )/f* = f D p w(x)D p e* kx dx = 

p=0 Jo 

= V^(— l) p / w(x)D 2p e tkx dx 

p=o Jo 

n •'O 


and together with (1.1.20) we recover the usual Fourier expansion we had before, namely 

1 C 


(1.1.23) 


w s (k) = w(k) = ^J w (0 e ik(d C 
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The completion of {e ,fcl } in ff J [0,2x) gives us the Parseval equality, compare (1.1.15) 


(1.1.24) 


\\w - svhIh* = E K(*0l 2 ll e, *1lH* = E [l^)l 2 - 27r E A;2p ] ^ 

\k\>N |fc|>^ P -0 

E^ 2p - 2 tt- E K^)l 2 = E^ 2p -ll^-^HI 2 - 

p= 0 \k\>N P= 0 


Since 


(1.1.25) 


Const 1 (l + jV 2 )* /2 < 



< Const 2 (l + IV 2 )?, 


we conclude from (1.1.24), that for any u;e/f'[0,27r) we have 
(1.1.26) \\w - SjvHI < Const, • ti>eff*[0,27r). 

Note that Const, = Consti ■ ||n; - S N w\\ H - — * 0. This kind of estimate is usually refered 

N—*o° m 

to by saying that the Fourier expansion has spectral accuracy , i.e., the error tends to zero 
faster than any fixed power of N, and is restricted only by the global smoothness of w(x). 

We note that as before, this kind of behavior is linked directly to the spectral decay of 
the Fourier coefficients in this case. Indeed, by Cauchy- Schwartz inequality 


(1.1.27) 


i -mi i - mi IMItf* • ll etfcx lltf* ^ 

K*)l = K(*)l < — < 


pifcx I 2 

e I H' 


(2*E‘ 0^)’ 




< Const 


(1 + W 2 )i 


In fact more is true. By Parseval equality 


Ml 


2 

H a 


OO OO / a \ 

E W*)l J lk“’llt. = 2x E E* Jp l*WI ! . 

ic~ — OO k=-oo \p =0 / 


and hence by the Reimann-Lebesgue lemma, the product (1 + |^| 2 ) ? |^(^)| is only bounded 
(as asserted in (1.1.27), but in fact it tends to zero, 


(i + |fc| 2 )*Hfc)|— ■* 0- 

fc— ►OO 

Thus, w(k) tends to zero faster than \k\-‘ for all w(x)eH a . This yields spectral convergence, 
for j i 

M - $Ml a = 27r E K fc )l 2 < Const. E n 4. Ifcpv - Copst - twtt 

|*:| > JV \k\>N t 1 + \ K \ ) 
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i.e., we get slightly less than (1.1.26), 


w — .SV^II < Const.- 


0 3 > 1. 


N‘~? N—*oo 

Moreover, there is a rapid convergence for derivatives as well. Indeed, if u>(£)ei/*[0, 27 t) 
then for 0 < <j < s we have 


= E ( 2 *'E k 2 p )\™W \ 2 

|*|>JV p= 0 

< Const. E ( X + l*| 2 rK*)| 2 
l*l>JV 


^ fritv-t V (^ 7r Sp=0 ^ 2P ) | . / , n , 2 

< Const. ^ — - \w{k)\ 


(! + N 2 Y- 


Hence 


< 


Const. 


|tu — ■SVHI/f* 
N 2 (‘~ v ) 


(1.1.28) (('to — 5Vu>||/r® < Const, • jj— — , a < s, w€H s [0,2tt) 

with Const, ~ ||tfl — SVHIh* * 0. Thus, for each derivative we “lose” one order in the 

N—*oo 

convergence rate. 

As a corollary we also get uniform convergence of S N w(x ) for #*[0, 27r)-functions u>(x), 
with the help of Sobolev-type estimate 


(1.1.29) Max 0 <x< 2 W |w( 2 :)| < Const.||v||tfi. 

(Proof: Write v(x ) = u(x 0 ) + J‘ 0 v'(x)dx with v(x 0 ) = ± / 0 2ir v(x)dx, and use Cauchy- 
Schwartz to upper bound the two integrals on the right.) 

Utilizing (1.1.20) with v(x) = w(x) — S^w(x) we find 

Maxo<x^ 2 ir l^(^) ^ f ^^(^){ ^ Const. ||iw 

(1.1.30) 

- Const ' aFTTTE °> 

JV N—>oo 


Corollary: Assume w(x)eH 1 [0 y 2 tt). Then 

oo 

(1.1.31) ro(x) = E w(k)e ,kx . 

k= — oo 
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In closing this section, we note that the spectral-Fourier projection, Snw(x), can be 
rewritten in the form 


(1.1.32a) 


S N w(x) = £ w(Jc)e ikx =—( *w({) £ e* x -Vd{ = 

k=-N J t=° k=-N 

= r d N (x - t)w(t)dt 

J(= 0 


where 

(1.1.324) D N (x -()=-{: e'«~» = ("+*)(»-« 

k=-N 27T s i n ^£^1 j 

Thus, the spectral projection is given by a convolution with the so-called Dirichlet kernel, 

(1.1.33) 

Now (1.1.23) reads 

(1.1.34) \w(x) - D n (x) * io(:e)| < Const, • j-, Const, ~ 


1 sin (N + 1) x 
D„ ( x)=- 


* X 

sin 2 
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1.2. The Pseudospectral (Collocation) Approximation 

We have seen that given the “moments” 

w(k) = C" w(()e- iki d(, —N < k < N, 

27 r o 

we can recover smooth functions w(x ) within spectral accuracy. Now, suppose we are given 
discrete data of u;(x): specifically, assume w[x) is known at equidistant collocation points 

(1.2.1) w v — w(x v ), x v = r-\-vh, v = 0, 1, • • • , 2N. 


Without loss of generality we can assume that r-which measures a fixed shift from the origin, 


satisfies 


( 1 . 2 . 2 ) 


0 < r < h = 


2 N + l 


Given the equidistant values w U) we can approximate the above “moments,” w(k), by the 
trapezoidal rule 2 


(1.2.3) 


h 2N+1 

<K k ) = r E " w - 




Using w(k) instead of w(k) in (1.1.7), we consider now the pseudospectral approximation 


(1.2.4) 


i> N w= X) 


The error, tu(x) — consists of two parts: 

w(x) - i> N w(x) = Y1 ™( k )e ikx + - w{k))e ikx . 

|fc|>N |fe|<// 

The first contribution on the right is the truncation error 


(1.2.5) 


Tnw(x) = (7 — Sn)'w(x) = EZ w(k)e ,kx . 

|fc|>^ 


We have seen that it is spectrally small provided w(x) is sufficiently smooth. The second 
contribution on the right is the aliasing error 


( 1 . 2 . 6 ) 


Anw(x) = ["&(&) — ih(fc)]e* fcx - 

\k\<N 


This is pure discretization error; to estimate its size we need the 

and ^ 2 " indicate summation with ^ of the first and respectively, the first and the last last terms. 
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Poisson’s Summation Formula (Aliasing) 


Assume Then we have 

(1.2.7) w(k) = Y e ip( ~ 2N ^ T w{k + p(2N + 1)). 


p— — oo 


Proof: For w(x)eH l [0,2‘K) we insert its Fourier expansion 


( 1 . 2 . 8 ) 


w(k ) = 


1 


2N 


Y^w(x v )e 


-xkx v 


1 


2N 


1 . f\ 


X) w(j)e ijx 


J~ oo 


" i / cXj / 


2iV + l^o ' ' 21V + ._ 0 

Since w(x) is assumed to be in H 1 i the summation on the right is absolutely convergent 
(1-2.9) Y < fe(l +i 2 )l^0‘)| 2 • Y 7T1-) < Const.|H| H i, 

i=-oo V j j 1 ' 7 ) 

and hence we can interchange the order of summation 

1 oo 22V 

(1.2.10) w(k) = T, *0) E 


»/=0 


Straightforward calculation yields 

j 2 N 


2 N + 


is 

n 




27V 


(1.2.11) 


— p i(}-k)r _ 


i/=0 


1 


21V + 


I-E.* 

1 i /=0 




377+T — 


r _ i 

2 - ^ +1 = 0 j — k ^ 0[mod 21V + 1] 


21V + 1 




l 21V + 1, 

and we end up with the asserted equality 


j - k =p-(2N + l). 


1 


2N 


™( k ) = J2 ™ti) • 2iV + - E) e ,(i fc)jv = E + p( 2N + !)) • e ,p ( 2W+1 ) r . 


J=-oo 




p= — oo 


We note that once ru(x) is assumed to be smooth, it is completely determined (pointwise, 
that is) by its Fourier coefficients w(k)] so are its equidistant values w„ = w(x u ) and so are 
its discrete Fourier coefficients w(k). The last formula shows that w(k) are determined in 
terms of w(k), by folding back high modes on the lowest ones, due to the discrete resolution 
of the moments of w(x ): all modes = &[mod21V + 1] are aliased at the same place since they 
are equal on the gridpoints 


( 1 . 2 . 12 ) 


0 i(k+p(2N+-i))x„ _ e ip(2N+l)r _ p ikx„ 
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Let us rewrite (1.2.7) in the form 

w(k) = w(k) + Y e’ p(2Ar+1)r • w(k + p(2N + 1)). 

p^o 

Returning to the aliasing error in (1.26), we now have 


(1.2.13) 


A n w(x ) = E 

\k\<N 


E e? p < 2N+1)r • w(k + p • (2 N + 1)) 
Lp^O 


Akx 


We note that T^w(x) lies outside 7 Tjv while Aatu;(x) lies in 7 Tjv, hence by 7P-orthogonality 
||xtr(=c) - xp N w(x)\\ 2 H . = E (1 + l*| 2 )' ' l™(*OI 2 truncation 

|fc|>AT 

+ Y, (1 + |A;| 2 )* • | E e tp ^ 2N+1 ^ r ■ w(k + p(2N + 1))| 2 <— aliasing. 

|fc|<jv p/o 

Both contributions involve only the high amplitudes - higher than N in absolute value; in 
fact they involve precisely all of these high amplitudes. This leads us to aliasing estimate 

Y (1 + l*| 2 )*l E e' p ( 2N+l > ■ w(k + p(2N + 1))| 2 < 


\k\<N 


p^O 


E E(1 + I* + P(2iV + 1)| 2 )X* + J>(2W + 1))|’ 


(1.2.15) 


•Max| fc |<^E 


i + i*r 


5SLi + l* + p-(2Jv + i)|*j 

1 + JV a 


< 


\\Tnw(x)\\ 2 h ..Y 

p^o L 


1+4 p 2 N 2 


Hence, we conclude that for any s > \ we have 
(1.2.16) 

Augmenting this with our previous estimates we end up with spectral accuracy as before, 
namely 

(1.2.17) ||u> - < Const, • weH s [0,27t), s > a > 1 

We observe that 'tp^w(x) is nothing but the trigonometric interpolant of u/(x) at the equidis- 
tant points x = x^: 


||i4 w i(;(x)||jf. < Const, • ||7Vu;(x)||tf., s > -. 


N 


tf} N w(x)\ x=Xii = Y 


k— — N L 


2N 


2 N + 


1 i /=0 




(1.2.18) 


2N 


N 


= Y 


v-0 


2 N + 


T E 

I 1 


0 ik{ii-v)h 


W 


M- 


k=-N 
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This shows that ip N is in fact a ^seudospectral projection, which in the usual sin-cos formu- 
lation reads 


(1.2.19) 


N 

'IpN'W = Y cos kx + bk sin kx 

k~0 



2 

2N + 1 


2 N 


Y w M 

u-o 


cos kx u 
sin kx u 


Thus, trigonometric interpolation provides us with an excellent vehicle to perform approxi- 
mate discretizations with high (= spectral) accuracy, of differential and integral operations. 
These can be easily carried out in the Fourier space where the exponents serve as eigenfunc- 
tion. For example, suppose we are given the equidistant gridvalues, w u , of an underlying 
smooth (i.e., also periodic!) function w(x) y w{x)tH*[ 0,27r). A second-order accurate discrete 
derivative is provided by center differencing 




Note that the error in this case is, 0(h 2 ) = u>( 3 )(£)/i 2 , no matter how smooth w(x) is. 
Similarly, fourth order approximation is given via Richardson procedure by 


dw 

dx 


(x 




8[ttJ„+i - UJ„-l] - [w v + 2 - W 1 /- 2 ] 

12 h 


+ 0(h A ). 


The pseudospectral approximation gives us an alternative procedure: construct the trigono- 
metric interpolant 


N 


( 1 . 2 . 20 ) 


'ipNw(x) = ^2 w(k)e 
k=-N 


ikx 


W 


(k) 


2N + 


1 2 N 

1 i /=0 


Differentiate - in the Fourier space this amounts to simple multiplication since the exponen- 
tials are eigenfunctions of differentiation, 

d N 

(1.2.21) —ipNw(x) = Y w(k)ike tkx , 

k=-N 


and we approximate 

(1.2.22) ~ X *') = ’^ } N'U){x)\x=x v + spectrally small error. 

Indeed, by our estimates we have for w(x)eH a [0 ) 2tt), s > 1 , 


(1.2.23) 


hlaXg<x<2x 



d_ 

dx 


ipNw(x)\ < Const. ||tn(x) — ipNw(x 


h 3 < 


Const, 

N- 2 
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which verifies the asserted spectral accuracy. Similar estimates are valid for higher deriva- 
tives. To carry out the above recipe, one proceeds as follows: starting with the vector of 
gridvalues, w = (w 0 , • • • , w 2 n), one computes the discrete Fourier coefficients 


(1.2.24) 


2N 


w(k ) = V* w v 

1 ' 2N + 1 to 


e ~ ikXv , —N < k < N, 


or, in matrix formulation 

(1.2.25) 

then we differentiate 

(1.2.26) 

or in matrix formulation 
(1.2.27) 


' w{-N) ' 


W 0 


= F 

\ 

. . 


l W 2N _ 


> F kv = 


1 


2N + 1 


0 -ikx v , 


w(k) — > ikw(k), 


' w(-N ) ' 


' w(-N) ' 


' -iN 


A 


, A = 

■ . 

w(N) 


w(N) 


iN _ 


and finally, we return to the “physical” space, calculating 


(1.2.28) 

or in matrix formulation 


N 


ikw(k)e lkXu , v = 0, 1, • ■ • , 2N, 

k=—N 



' £<*•> ' 


- -iNw(-N) ' 

(1.2.29) 


= F* ■ (2N + 1) 

; 


. . 


iNw(N) 


, (2N + l)FZ k = e' kx ''. 


The summary of these three steps is 
(1.2.30) 


w'(x 0 ) 


Wo 


= ipD 

: 

_ w\x 2N ) _ 


. W2N . 


, ipD = (2N + 1)F*AF, 


where xpD represents the discrete differentiation matrix, and similarly ipD” for higher deriva- 
tives. 

Note : Since (2 N + 1 )F*F = I 2 n+i (interpolation!) we apply 'tpD 3 = (2N + 1)F*A‘ , F’. How 
does this compare with finite differences and finite-element type differencing? 
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In periodic second-order differencing we have 


fourth order differencing yields 


0 1 

-1 0 


FD 2 = — 


8 -1 


-8 0 


0 -1 

0 

0 1 

-1 0 


1 -8 

1 


-1 

8 -1 


-1 

0 8 

1 -8 0 


In both cases the second and fourth order differencing takes place in the physical space. 
The corresponding differencing matrices have finite bandwidth and this reflects the fact that 
these differencing methods are local . Similarly, finite-element differencing, 


e^-i + 6 wl + 


corresponds to a differencing matrix 


Wy+l ~ Wy-l 

2h 


0 1 

-1 0 


-1 0 


We still operate in physical space with O(N) operations (tridiagonal solver) and locality 
is reflected by a very rapid (exponential decay) away from main diagonal. Nevertheless, if 
we increase the periodic center differences stencil to its limit then we end up with global 
pseudospectral differentiation 


(1.2.31) 


s JjjT+rE’ 


kx u \ ikx v 


recall the Dirichlet kernel (1.1.33) 


ikx — iNx 

e — e 


gt( 2 Ar+i)x _ j sin(JV + i)a 
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and its derivative, 


y- d sin(N + \)x _ (N + l)cos(JV+ £)xsin§ - |cosfsin(AT + |)i 

(1.2.33) ^,ke = - — - S? | 


T. ike 

k=-N 


so that 

(1.2.34) 

Hence (1.2.31), (1.2.34) give us 

(1.2.35) 


*_»)» _ (w + |)c«[(iy + !)(■/ -<«)t] 




w 


A 2N 1 

'K) S d^ NW ^ = S o-7^V • h^U = 


-2sin(^) 


2sin(^^) 

In this case ifiD is a full (2JV + 1) x (2 AT + 1) matrix whose multiplication requires 0(N 2 ) 
operations; however, we can multiply ^ D[w ] efficiently using its spectral representation from 
(1.2.30), 

i)D = (2N + l)F*AF. 

Multiplication by F and F* can be carried out by FFT which requires only 0(N log N) 
operating and hence the total cost here is almost as good as standard “local” methods, and 
in addition we maintain spectral accuracy. 

We have seen how the pseudospectral differentiation works in the physical space. Next, 
let’s examine how the standard finite-difference/ element differencing methods operate in 
the Fourier space. Again, the essential ingredient is that exponentials play the role of 
eigenfunctions for this type of differencing. To see this, consider for example the usual 
second order center differencing, D 2 {h) y for which we have 


(1.2.36) 


D 2 {h)e 


ikx | 


e ik Xv+ i _ e isin(Jfch) 


,ikx | 


c— Xj, 2^ fa (*=*., 

The term — - ?C - £ 2 is called the “symbol” of center differencing. By superposition we obtain 
for arbitrary grid function (represented here by its trigonometric interpolant) 


N 


(1.2.37) 

that 

(1.2.38) 


ipNw(x) = w(k)e' 

k=-N 


ikx 


w v +i - iy. i 

2 h 


N 


= D 2 (h)ip N w = y w(k)D 2 {h)e 

k=-N 


ikx\ 


| x—x u 


N 

= E 

k=-N 


isin(kh) ik 

^ L w(k)e tlcx \ 
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It is second-order accurate differencing since its symbol satisfies 


(1.2.39) 


i s'm(kh) 


ik + 0(k 3 h 2 ). 


Note that for the low modes we have 0{h 2 ) error (the less significant high modes are differ- 
enced with 0( 1) error but their amplitudes tend rapidly to zero). Thus we have 


(1.2.40) 


|| ±4,„ w - DAWvwf = E I* _ 

ax k=-N n 

< Const. h A ^ (1 + |ifc| 2 ) 3 |^)| 2 < Const. /t 4 • HV'tf'HIff 3 > 
\k\<N 


and this estimate should be compared with the usual 

— w{xj) — Wl/ -~ < Const. h 2 • Max x |u/ 3 )(x)| . 

dx 2h 11 

The main difference between these two estiamtes lies in the fact that (1.2.41) is local, i.e., 
we need the smoothness of w(x) only in the neighborhood of x = x v and not in the whole 
interval. The analogue localization in the Fourier space will be dealt later. 

Similarly, we have for fourth order differencing the symbol 


.1 I" sin** sin2kh 

i— 4 

3 h 2 h 


= ik + 0(k*h 4 ). 


In general, we have difference operators whose matrix representation, D, is of the form 


(1.2.41) 


D = [d jk ] -N<j,k<N 


and it is periodic and antisymmetric (here [£] = £[mod 2 N + 1]) 


[1.2.42) 


(i) periodicity dj k = d[ k _j] 

(ii) antisymmetry d jk = — d k j . 


Matrices satisfying the periodicity property are called circulant, and they all can be diago- 
nalized by the unitary Fourier matrix 


(1.2.43) 


D = U*AU, U = (2N + l)i - F, U*U = I 2N + 
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Indeed, with p — q = i we have 


W’DUU = 2?rp[ E 1 = 


Z[p-g]* 


p,q=-N 
N 


= TfTTT £ e«'+<«+'Wd M e-“('+-‘) 

ZJV 1 l,q=-N 


(1.2.44) 


1 i E e«% • E 

1 # — * r q=-N 


2N + 1 a?.* 

( o j^k, 


+g*) 


= < 


N 


E i = *> 

l=-N 


and using the antisymmetry we end up with symbols A* 


(1.2.45) 


N 

A = diag(A_//, • • • , A m), A* = 2i E sin(WA). 

/=! 


As an example, we obtain for the finite-element differencing system 


.s'mkh /4 1 ■ 


(1.2.46) 


( 


6i 


6 

sin(fc/i) 


+ -e“" + 


h 4 + 2 cos(fc/i) 


r-“) = 

= ik + 0(A 4 ). 


In general, the symbols are trigonometric polynomials or rational functions in the “dual 
variable,” kh, which has “exact” representation on the grid in terms of translation operator 
(polynomials or rational functions), and accuracy is determined by the ability to approximate 
the exact differentiation symbol ik for | A: | ~ 1, see Figure 2. 
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We have seen that the spectral and the pseudospectral approximations enjoy what we 
called “spectral accuracy” - that is, the convergence rate is restricted solely by the global 
smoothness of the data. The statement about “infinite” order of accuracy for C°° functions 
is an asymptotic statement. Here we show that in the analytic case the error decay rate is 
in fact exponential . 

To this end, assume that 

(1.3.1) w(z ) = ^2 w(k)e lkx , \Im z\ < rj < r) 0 , 

k— — oo 

is 27r-periodic analytic in the strip —770 < Im z < rj 0 . The error decay rate in both the 
spectral and pseudospectral cases is determined by the decay rate of the Fourier coefficients 
w(fc). Making the change of variables ( = e xz we have for 

(1.3.2) u(() — w(z — +i£n(), 
the power series expansion 

00 

(1-3.3) v(() = £ *(4)C‘. 

k= — OO 
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By the periodic analyticity of w(z ) in the strip \Imz\ < t) < rj 0 ,v(Q is found to be single- 
valued analytic in the corresponding annulus 

(1.3.4) e -T?0 < \rj\ < e tJO , 
whose Laurent expansion is given in (1.3.3): 

(1.3.5) w(k) = - j v(C) C (k+1) dC, e"” 0 <r<e'*. 

' " Z7TZ J j£j=r 

This yields exponential decay of the Fourier coefficients 

(1.3.6) |tu(fc)| ^ M(Tj)e~ kv , M(rj) = Max i/mi |< T> |iw(z)| > 0 < rj < rj 0 . 

We note that the inverse implication is also true; namely an exponential decay like (1.3.6) 
implies the analyticity of w(z). Inserting this into (1.1.24) yields 

II 1X7 


(1.3.7) 


S N w\\ 2 = 2tt 2 Hk)| 2 < 

\k\>N 

< 2tt • M 2 { rj) • £ e~ 2hn = 2.^1 • 




e 2lJ — 1 


and similarily for the pseudospectral approximation 


(1.3.8) 


w — ip N w || 2 < Const. 


M 2 (v ) 


Nr, 


e 2fl — 1 

Note that in either case the exponential factor depends on the distance of the singularity 
(lack of analyticity) from the real line. For higher derivatives we likewise obtain 

,—2Nt) 


(1.3.9) 


w — 5jvttf||jf<r -f 1 1 — tPn w \\‘h<’ — Const. N 2<r ■ M 2 (r]) 


e 2n _ i 


We can do even better, taking into account higher derivatives, e.g., 

hL d £w- ld( - 


(1.3.10) 

so that with 

(1.3.11) 

we have 

(1.3.12) 

and hence 

(1.3.13) 


= e 2 *^Max| C | =e ,|u (j) (OI, 

3=0 


fc|w(Jb)| < M 1 (T ] )e 


-kt) 


,-2 Nr, 


| to — Snw\\ 2 h „ + ||tw — ip N w\\ 2 Ha < Const .Ml(rj) — — - 
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1.4. The Non-Periodic Problem - Chebyshev Approximation 


We start by considering the second order SL problem 

(1.4.1) — n/I — x 2 — (\/l — x 2 — ^V’( I )> — 1 < x < 1. 

This is a special case of the general SL problem 

(1.4.2) XV’ = 7-rj- (p( x )t~') + = W( x )> P,<l, w > 0. 

(V{x) ax \ ax j uj[x ) 

Noting the Green identity 

(1.4.3) o,( x ) = / -(pV0V + = P{ x )[i>^]\ b a + (V’» -^^)u.(i)> [V’, = V’V’' - V'V’'. 

we find that L is formally self-adjoint provided certain auxiliary conditions are satisfied. In 
the nonsingular case where p(a) • p(b) ^ 0, we augment (1.4.2) with homogeneous boundary 
conditions, 


(1.4.4) 'ip(a) ~ <f>(a) — 0, ^(b) = <j)(b) = 0. 

Then L is self-adjoint in this case with a complete eigensystem (A*, ^(x)): 
tu(as)eL w ( x )[a, &] has the “generalized” Fourier expansion 

( 1 . 4 . 5 ) w(x) ~ x. n k ) = — — 


each 


1 f h 

w( k ) = p-p- J * w(x)^ k (x)uj(x)dx. 


with Fourier coefficients 

( 1 . 4 . 6 ) 

The decay rate of the coefficients is algebraic: indeed 

w(k) = 


(1.4.7) 


1 

l 

IML 

' y k { 

i 

i 

Ml 

' a7 ( 

i 

II./. 1 1 0 

(p(x) 


(p(x) • (V’i.HlS + (V'fc, 


J-l 


X V W|1 LTM — / I a 1 \ J \ r 7 — y 

U'¥'*||« j=0 A A: "'it 

The asymptotic behavior of the eigenvalues for nonsingular SL problem is 

I 2 

7rk 




Const, k 2 
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and hence, unless w(x ) satisfies infinite set of boundary restrictions we end with algebraic 
decay of w(k ) 

1 K*) ././ i& Const * 


w(k) 


HV'fcl 


V’UxMx)!, 


™ v ' lo k 2 

This leads to algebraic convergence of the corresponding spectral and pseudospectral pro- 
jections. 

In contrast, the singular case is characterized by having, p(a) = p(6) = 0, in which case 
L is self-adjoint independently of the boundary conditions since the brackets [ , ] drop, and 
we have spectral decay-compare (1.1.22) 


(1.4.8) 


M.k) = — i— - • A. . (*,£(•)«,)„ < 

HAUL K HAL 


provided w(x) is smooth enough; that is, the decay is as rapid as the smoothness of tt;(x) 
may permit. 

Returning to the singular SL problem (1-4.1) we use the transformation 

1 d 


(1.4.9) 

which yields 

(1.4.10) 


d Id 

X = COS 9, "jT" = ~dx "j/T 


dx g do y/1 - X 2 dd 


d 2 


dd 2 


(f)(9) = A (f)(9), (f)(9) = V>(cos 9), 


obtaining the two sets of eigensystems 

(1.4.11a) (X k — k 2 ,(f> k = cos k9), 

and 

(1.4.116) (A* = A: 2 , (f> k = sin k9). 

The second set violates the boundedness requirement which we now impose 


(1.4.12) |V>*(±1)| ^ Const., 
and so we are left with 

(1.4.13) (A* = k 2 ^ k (x) - cos(kcos~ 1 x)). 
The trigonometric identity 

cos(A: + 1)9 — 2 cos 9 cos k9 — cos (k — 1)0 
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yields the recurrence relation 


(1.4.14) i>k+i{x) = 2xip k (x) - V'fc-ifa), V'o(x) = 1 , V’i( a: ) = x > 
hence, V’fc( a: ) are polynomials of degree k - these are the Chebyshev polynomials 

(1.4.15) Tfc(x) = cos (A: cos -1 x) 


n . in , TMTf \ \ [' T k (x)Tj{x) 

(1.4.16) (?*(*), Tj{x)) w = ^i_ x r dx 


= (l-x 2 )-*, 


' 0 

j ^ 

< 1 \TkWl = \ 

j = k > 0, 

. Iirolll = ^ 

j = k = 0. 


In analogy with what we had done before, we consider now the Chebyshev-Fourier expansion 

(w(x),T k (x)) u 


II T k 


(1.4.17) w(x) ~ ™(k)Tk(x), w(k) 

k= 0 

To get rid of the factor j for k = 0 we may also write this as 

oo 

(1.4.18a) w(x) ~ 

A :=0 


where 

(1.4.186) 


w(Jc) = 


(w(x) } T k (x)) u , 

7r/2 


2 f 1 w(x)cos(kcos 1 x)dx 

7T 7-1 \/l ~ X 2 


or using the above Chebyshev transformation 


(1.4.18c) 


2 f* 

w(k ) = — / u;(cos £) cos 

7T J£=0 


Thus, we go from the interval [—1,1] into the 27r-periodic circle by even extension, with 
Fourier expansion of u;(cos0), compare (1.1.9) and see Figure 3, 

1 /* 27r 2 f* 

w(k ) = — / u;(cos £) cos = — / u/(cos £) cos k(d£. 

7T 7^0 7T */£=0 

Another way of writing this employs a symmetric doubly infinite Fourier-like summation, 
where 


(1.4.19) 


1 00 

™( x ) ~ 9 £ ™(k)T k (x) 

Z h — — oo 
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w(9) = w(x — cosO ) 



with T_t(x) = Tk(x ) and 

(1.4.20) w(k) = - f 1 W ^ T t&dx, -oo < k < oc 

7T J-l Vl — X 2 

The Parseval identity reflects the completeness of this system 

(1.4.21) 

which yields the error estimate 

Ik - S N w\ \\. = k(*)| 2 - 

* k>N 

Now, in order to get a measure of the spectral convergence we have to estimate the decay 
rate of Chebyshev coefficients in terms of the smoothness of w(x) and its derivatives; to 
this end we need Sobolev like norms. Unlike the Fourier case, {T*.(x)} is not complete with 
respect to H a - orthogonality is lost because of the Chebyshev weight. So we can proceed 
formally as before, see (1.1.24), 

(1.4.22) |k - Sswfr = 2tt £ |t&(*)| a < £ 77 t ^ kWI* 

k>N k>N V A t" iv ) 

i.e., if we define the Chebyshev- Sobolev norm 

iMik = £(i + 

k= 0 


*ko)i>+i5:k*)i’ 

Z Jfc/0 


= jEW*)P 

% fc= o 
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then we have spectral accuracy 


1 1 to — Snw\\t < Const, • — - , weH^[— 1,1], 

N 3 

In fact the H? space can be derived from appropriate inner product in the real space as done 
in Fourier expansion. The correct inner product is given by - compare (1.1.19) 

* 3 r 2x d? p d? p 

(1.4.23) (w 1 ,w 2 ) H 2.=Y1( LPw '> LPw 2)t^Y;J ^^Wi(cos0 )-— iu 2 (cos e)d6, 

p=° x-cosep =0 e ~° 


(1.4.24) 


0 j / k, 

{Tk, Tj) H 2. = < 7T 3 

— 2J & > j = k (with 7r factor at j = k = 0). 
. 2 p=o 


Hence the Fourier coefficients in this Hilbert space behave like 


(1.4.25) 


(u;(x), T fc ) H a. ~ 53(1 + k 2 ) 2, w{k), 


and the corresponding norm is equivalent to 


(1.4.26) 


Hlk* ~ TXi + ^ 2 ) 2 *l^(^)l 2 - 


The reason for the squared factors here is due to the fact that L is a second order differential 
operator, unlike the first-order D = j- in the Fourier case, i.e., 


(1.4.27) 


E(i + W 2 ) 2 X*0l ! ~ E \\ lt ™\\t 


involves the first 2s-derivatives of w(x) - appropriately weighted by Chebyshev weight. This 
completes the analogy with the Fourier case, and enables us to estimate derivative as well- 
compare (1.1.28 - 1.1.29), 


(1.4.28) 


\w — Snw\\h* < Const,—— , a < s y weH 3 T [- 1,1]. 


Next, let’s discuss the discrete setup. Because we seek an even extension of the upper semi- 
circle we consider the case of even number of grid points - equally distributed along the unit 
circle. One choice is to look at 


e v = r + vh = (i/ + 1)-^- 


T 7 r h 

h = N' r= 2 
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Xu = COj[( V + j)£] 


Figure 4: 

which gives the usual Chebyshev interpolation points, see Figure 4. The second choice - 
which we will concentrate from now on - takes into account the boundaries, where we look 
at 

„ , 2tt . 

6 u = vh = v— (h=—,r = 0) 

which yields, see Figure 5, 

(1.4.29) x v = cos(uh), i/ = 0,l 

The trapezoidal rule is replaced by Gauss-Lobatto rule for (1.4.20), i.e., starting with 



x v = cos(uh) 
Figure 5: 


24 


(1.4.18c), we have 


(1.4.30) tt?(fc) = — / w( cos £) cos k(d( — > ► — Y] "w( cos 9 V ) cos k9 v * — , 

7T J f=0 7T iV 

and we end up with the discrete Chebyshev coefficients 


(1.4.31) 


>(*) = TT E VJ4 0<k<N. 

** is=0 


Indeed, this corresponds to the case of Fourier interpolation with even number of equidistant 
gridpoints 9 V in (A. 1.2), for 


2 N 

= 5- X) "' u v e ‘ 

^ n 




2 " 

= —J2"w l/ cos(kd u ). 

** u=0 

Then one may construct the Chebyshev interpolant at these N + 1 gridpoints 


(1.4.32) 


ip N w(x) = "w(k)T k (x). 


We have an identical aliasing relation 


(1.4.33) 


w(k) = ^ w(k + 2 pN) 


(compare A. 1.5), and hence spectral convergence, i.e., (compare (1.2.14) and (1.2.12)) 


(1.4.34) 


\w(x) - ip N w(x)\\ H * < Const, • — — , weHy, s > a, 


where Const, ~ ||iu||/f* . 


Example : We have the Sobolev imbedding of L°° space 




< Const ff • IMI/re, a > i. 


Consequently, 


Max x |u;(x) - Tp N w(x)\ < (Const, ~ ) • -77717, * > a > 1 
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In particular, with s - N + 1 we obtain the ” usual” estimate for the near minmax approxi- 
mation ( collocated at at x v = cos ^ 1 / + \)v)) 

e~ N 

Max x |m(x) — TpNw(x) \ < Const. ||u>|| H Ar+i • 

We briefly mention the exponential convergence in the analytic case. To this end we employ 
Bernstein’s regularity ellipse, E r , with foci ±1 and sum of its semi axis = r. Denoting 

(1.4.35) M(rj) = Ua.x xeEr \w(z)\, r = e n . 

We have 

Theorem. Assume io(x) is analytic in [-1,1] with regularity ellipse whose sum of semiaxis 
= tq = e'* ) > 1. Then 

||w(x) - 'ipNw(x) ||h* + _ < C° ns t. — • N 2tr e 2Nr) . 


Proof: The transformation 2 = (( + C _1 )/2 takes E TQ from the z- plane into the annulus 
Tq 1 < K| < tq in the £-plane. Hence, u(C) = 2 w [z - admits the power expansion 

(1.4.36) u(C) = 2 w ( ^y = f] w( k )( k , Tq 1 < Id < r 0 = e Co ; 

indeed, setting ( — e ie and recalling w(-k) = i b(k), the above expansion clearly describes 
the real interval [-1,1] 

OO 

(1.4.37) w(z = cos 8) — 'w(k) cos kO. 

k= 0 


By Laurent expansion in (1.4.36) 

( 14 - 38 ) - 2 hL„ e ~’ K<r< 

hence 

(1.4.39) m^)| < M{r})e ~ kj1 

and the result follows along the lines of (1.3. 7-8). 


Finally, we conclude with discussion on Chebyshev differencing. Starting with grid values 
w v at Chebyshev points x v — cos (vjf), one constructs the Chebyshev interpolant 

N 2 N 

( 1 . 4 . 40 ) "«<*) • r fc (x), w(k) = —J2"w u cos(kx„). 

k = 0 v=0 
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One can compute w(k), 0 <k<N, efficiently via the cos-FFT with 0(N log N) operations. 
Next, we differentiate in Chebyshev space 

( 1 . 4 . 41 ) ^-ip N w(x) = Y 2 

In this case, however, Tjt(x) is not an eigenfunction of instead ^Tk(x) - being a polyno- 
mial of degree < k — 1, can be expressed as linear combination of {Tj(x)}*Zo 0 n ^ act ^fe( x ) 
is even/odd for even/ odd fc ; .s): with Cq — 2, c^> o = 1 we obtain 

(1.4.42) S r ‘ (x)= l S T ' (X) ’ 

aj ' 0<j<fe 

fc-j odd 

and hence 

(1.4.43) ^-$ N w(x) = Y2 w(k)— T A X )- 

dx k _ 0 Ck o <j<k 

k-j odd 

Rearranging we get 

(1.4.44) ih'( fc ) = - £ P*(p) 

ox fc=0 C fc p >fc + i 

p+fc odd 

and similarly for the second derivative 

(1.4.45) w"(k) = — YL P(P 2 ~ ^ 2 )^(p)- 

C k P>k + 2 

p+J* even 

The amount of work to carry out the differentiation in this form is 0(N 2 ) operations which 
destroys the IVlogiV efficiency. Instead, we can employ the recursion relation which follows 
directly from (1.4.44) 

(1.4.46) w\k + 1) = w\k — 1) • Ck - 1 — 2kw(k). 

To see this in a different way we note that 

sin(A: + 1)8 — sin (k — 1)9 -f 2 sin 0 cos k8, 
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as asserted. In general we have 


(1.4.47) 


+ 1) = w { ’\k - 1) C4 _! - 2kw ( - a ~ 1) (k). 


With this w(k) can be evaluated using <9 (TV) operations, and the differentiated polynomial 
at the grid points is computed using another cos-FFT employing 0(N log TV) operations 


(1.4.48) 


d N 

tpN w {, x )\ x =x v = 53 "w'(k)cos kx v , 
k= 0 


dx 


with spectral/exponential error 
(1.4.49) 


Max x=Ii/ 1 ^^^(i) ^ ^ 


Const, • 


N’-<r 

Const, • e~ Nr) 


2 <,J<a ' 


The matrix representation of Chebyshev differentiation, D T , takes the almost antisymmetric 
form 

C)f. Xj Xfc 


{D T )jk = < 


-2(T% .n 

27V 2 + 1 


27V 2 + 1 


3 = k = 0 , 


. j = k = N. 
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2. TIME DEPENDENT PROBLEMS 


2.1. Initial Value Problems of Hyperbolic Type 


The wave equation, 

(2.1.1) w tt = a 2 w xx , 

is the prototype of P.D.E.’s of hyperbolic type. We study the pure initial- value problem 
associated with (2.1.1), augmented with 27r-periodic boundary conditions and subject to 
prescribed initial conditions, 

(2.1.2) w(a :,0) = /(x), w t (x,0) = g(x). 


We can solve this equation using the method of characteristics , which yields 


(2.1.3) 


f(x + at) + f(x - at) 
w(x,t) = + 


\ rx+at 

2 a J x—at 


g(s)ds. 


2 Zd Jx—at 

We shall study the manner in which the solution depends on the initial data. In this context 
the following features are of importance. 


1. Linearity: the principle of superposition holds. 

2. Finite speed of propagation: influence propagates with speed < a. This is the essential 
feature of hyperbolicity. In the wave equation it is reflected by the fact that the value of 
w at (x, i) is not influenced by initial values outside domain of dependence (x—at, x-f at) 

3. Existence for large enough set of addmissible initial data: arbitrary C initial data 
can be prescribed and the corresponding solution is C £°. 


4. Uniqueness: the solution is uniquely determined for — oo < t < co by its initial data. 

5. Conservation of Energy . The wave equation (2.1.1) describes the motion of a string 

with kinetic energy, J w 2 dx y and potential one, jT / w 2 x dx, (T / p = a 2 ). In order 

to show that the total energy 


E< 


Total 


= ?'/( 


2 2 2 
w, -f a w x 


)da 


is conserved in time we may proceed in one of two ways: either by the so called energy 
method or by the Fourier method. 
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The Energy Method. 

Rewrite (2.1.1) as a first order system 


(2.1.4a) 


d_ 

dt 


Ui 


' 0 

a 2 ' 

d 

Ui 


111 

_ u 2 


1 

0 

dx 

. U2 . 

) 

U 2 


dw 

~dt 

dw 

L dx - 1 


or equivalently, 

(2.1.46) 

The essential ingredient here is the existence of a positive symmetrizer, H > 0, 


du du 
dt = A fc' 


HA 


0 a 2 
a 2 0 


= A = A], H = 


1 0 
0 a 2 


(2.1.5) 

so that multiplication by H on the left gives 

(2.1.6) Hu t = A g u x . 

Now, multiplying by u T we are led to 

(2.1.7) (it, Hu t ) = (it, A,u x ), 

and the real part of both sides are in fact perfect derivatives, for by the symmetry of H , 


Re(tt, Hut) = l(u, Hut) + u ) = 

1 1 ^ 

~(u,Hu t )+ ~{ut,Hu ) = ^ 




and similarily, by the symmetry oM J3 we have 


11 & 
Re(it, A,U X ) = -(it, A,u x ) + -( A,u x , it) = 


(u,A s u) 


2\-> -*/ ■ -/ Qx [ 2 ' 

Hence, by integration over the 27r-period we end up with energy conservation, asserting 

t t n 

(2.1.8) — f (wl + a 2 wl)dx = — f (u, Hu)dx = f —{u,A s u)dx = 0. 

dt Jx dt Jx Jx ux 

We note that the positivity of // was not used in the proof and is assummed just for the 
sake of making (u, Hu) an admissible convex “energy norm.” 
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The Fourier Method. 


Fourier transform (2.1.4b) to get the O.D.E. 

dii 


njj-(M) = ikAu(k,t), 


(2.1.9) 

whose solution is 

(2.1.10) u(k,t) = e' kAt ii{k,0), 

where ii(k, 0) is the Fourier transform of the initial data. Now, for 

(2.1.U) 

we find 

( 2 . 1 . 12 ) 


A = TAT -1 , A 


— a 


—a a 

, T = 

1 1 

a 

? 


u l 


(k,t) = T u(k, 0) 


or 


(2.1.13) 


T~ x u{k,t) = 


, — ifca£ 


o 


o 


n ikat 


T _1 u(fc, 0) 


and hence (since the diagonal matrix on the right is clearly unitary), the L 2 -norm of 
T -1 u(A:, t) is conserved in time,i.e., 


(2.1.14) 


||T _1, u.(A:, t)|| 2 = ||r - 1 u(fc, 0)|| 2 , T-' = -~ 


1 —a 

— 1 — a 


Summing over all modes and using Parseval we end up with energy conservation 


(2.1.15a) 


[, 2 , 2 2xr a 2 f (w t -aw x y fw t + aw x y 

y H + a w,)dx = 4a J _ (— —) + dx 

= 4 a 1 J ||T-’u||Vi = i)f = Const. 


as asserted. 

We note that the only tool used in the energy method was the existence of a positive sym- 
metrizer for A, while the only tool used in the Fourier method was the real diagonalization 
of A\ in fact the two are related, for if A — TAT" 1 then for H = (r _1 )*T“ 1 > 0 we have 

(2.1.156) HA = (T~ 1 yAT~ l = A t = A?, A real diagonal. 

Energy conservation implies (in view of linearity) uniqueness, and serves as a basic tool 
to prove existence. It will be taken as definition of hyperbolicity . It implies and is implied 
by the qualitative properties (1) - (4). 
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We now turn to consider general P.D.E.’s of the form 
(2.1.16a) 


Qxl ^ 0 

-jfi = p ( x , t, D ) u , p i x , t,D) = J2 A i( x > 


J=1 

with 27r-periodic boundary conditions and subject to prescribed initial conditions, 

(2.1. 16£) u(x } 0) = /(z). 

We say that the system (2.1.16a) is hypberbolic if the following a priori energy estimate 
holds: 


(2.1.17) 


K*.0lk(«) ^ Const.T • |K*,0 )|| Mx)i —T <t<T. 


As we shall see later on this is equivalent to energy conservation with appropriate energy 
renorming. 

Here are the basic facts concerning such systems. 

The Constant Coefficients Case. 


Oil ^ Q 

(2.1.18) — = P(D)u, P(D ) = V Aj— , Aj = constant matrices. 

ut ox 


Define the Fourier symbol associated with P(D ): 


(2.1.19) 


P(ik) = i Ajkj, k = (k lt k 2 , • • ■ , k d )eR d , 

J-l 


which naturally arises as we Fourier transform (2.1.18), 

0 


—u(k, t) = P(ik)u(k,t). 


( 2 . 1 . 20 ) 

Sovling the O.D.E. (2.1.20) we find, as before, that hyperbolicity amounts to 

(2.2.21) ||e #(ifc)< || < Const T , —T < t < T, for all Jfc's. 

For this to be true the necessary Garding-Petrovski condition should hold, namely 

(2.1.22) |ReA[P(*Jfc)]| < Const. 

Example: For the wave equation, A[P(zfc)] = ±ika. 

But the Garding-Petrovski condition is not sufficient for the hyperbolic estimate (2.1.17) as 
told by the counterexample 

d_ 

dt 


«1 


a 

1 ' 

d 

Ui 

U2 _ 


0 

a 

dx 

u 2 
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As before, in this case we have A[P(z&)] = ±ika , hence the Garding-Petrovski condition is 
fulfilled. Yet, Fourier analysis shows that we need both ||^i(x, 0)||^ 2 ( x ) and || ^-(x, 0)||l 2 ( x ) m 
order to upperbound ||?x 1 (x, i)||i, 2 ( x ). Thus, the best we can hope for with this counterexample 
is an a priori estmate of the form 

IKx, OIIm*) - Const T ■ ||u(x, 0)||//i( x ), -T <t <T. 

We note that in this case we have a ”loss” of one derivative. This brings us to the notion of 
weak hyperbolicity. 

We say that the system (2.1.16a) is weakly hyperbolic if there exists an 5 > 0 such that the 
following a priori estimate holds: 

IKMJItaw ^ Constr • ||u(x,0)||h*(x), -T <t<T. 

The Garding-Petrovski condition is necessary and sufficient for the system (2.1.18) to be 
weakly hyperbolic. The necessary and sufficient characterization of hyperbolic systems is 
provided by the Kreiss matrix theorem. It states that (2.1.21) holds iff there exists a positive 
symmetrizer H(k ) such that 

(2.1.23) R e[H(k)P(ik)} = 0, 0 < m < H(k) < M, 

and this yields the conservation for ||u(x, f)|| = 

27cJ2{u(Kt)),H(k)u(k,t)) 

k 

is conserved in time. 

Remark : For an a priori estimate forward in time (0 < t < T), it will suffice to have 

(2.1.24) Re[tf(ifc)P(iJfc)] = hH{k)P(ik) + P{ik)H(k)} < 0. 

Indeed, we have in this case 

~(u(k),H(k)u(k)) < (R e lH(k)P(ik)]u(k),u(k))<0, 

and hence summing over all k’s and using ParsevaPs 

HM)Hi 2 (x) < ~ ll«(*.0)ll£ a (,)- 

Special important cases are the strictly hyperbolic systems where P(ik ) has distinct real 
eigenvalues, so that P(ik ) can be real diagonalized 

P(k ) = *r(ifc)ACJb)r“ 1 (Jb), 
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and, as before, H(k) = (T -1 ^))*!’ -1 ^) will do. The other important case consists of 
symmetric hyperbolic systems which can be symmetrized in the physical space, i.e. there 
exists an H > 0 such that 

HAj = A it = Aj a . 

Most of the physically relevant systems fall into these categories. 


Example: Shallow water equations (linearized) 


d_ 

dt 


u 

d 

U 

t d 

u 

V 

+ A i — 

V 

+ a 2 ^~ 

V 

. 4 > . 

dx 

.t. 

dy 

A . 


= 0, 


with 


A\ = 

can be symmetrized with 


U 0 

0 

1 


Vo 

0 

0 

0 

u 0 

0 

j ^2 — 

0 

Vo 

1 

. <t> 

0 

u 0 


0 

</>0 

Vo . 


H 


4>o 


<t>a 


The Variable Coefficient Case. 


(2.1.25) 


~ = P(x, t, D)u. 


This is the motivation for introducing the notion of hyperbolicity as is: freeze the coefficients 
and assume hyperbolicity of u t = P(x 0 ,t 0 ,D)u uniformly for each (x 0 ,t 0 ) ; then (unlike the 
case of weak hyperbolicity), the variable coefficients problem is also hyperbolic. 


Remark : This result is based in the invariance of the notion of hyperbolicity under low- 
under perturbations; it restricts hyperbolic system to be of first-order. 

So far we have dealt with hyperbolicity via the Fourier method studying the algebraic 
properties of its symbol P{ik)\ we can also work with the energy method . 


For example, if we assume that P(x,t, D ) is semi-bounded , i.e., if 
(2.1.26) -M\\u\\ 2 Lt{x) < Re(u, P{x, t, D)u) L , {x) < M||u||^ (l) , 0 < M, 

then we have hyperbolicity. 
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Example: The symmetric hyperbolic case Aj(x y t) = we can rewrite such sym- 

metric problems in the equivalent form 

du 1 , du / A \ , T> T~> ^ 

= H + ' b = -2?&7' 

In this case the symmetry of the Aj’s implies that \ Yjj A/Jjj + Sj skew- 

adjoint, i.e., integration by parts gives 




Therefore we have 

Re(tir, P^Xy tj -0)^)l 2(®) = Re(J?iij ^)l 2 (x)> 

and hence the semi-boundedness requirement (2.1.26) holds with M = ||ReB||. Conse- 
quently, if Aj(x } t) are symmetric (or at least symmetrizable) then the system (2.1.16a) is 
hyperbolic. 
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2.2. Initial Value Problems of Parabolic Tipe. 


The heat equation, 

(2-2.1) u t — au xx , a > 0, 

is the prototype of P.D.E.’s of parabolic type. We study the pure initial-value problem 
associated with (2.2.1), augmented with 27r-periodic boundary conditions and subject to 
initial conditions 

(2-2.2) u(x,0) = /(x). 

We can solve this equation using the Fourier method which gives 
(2-2.3) u(M) = e~ akH f(k). 

It shows the dissipation effect (= the rapid decay of the amplitudes, \u{k, t)|, as functions of 
the high wavenumbers, |fc| > 1) in this case, which is the essential feature of parabolicity. 
As before, we study the manner in which the solution depends on its initial data. 

1. Linearity : the principal of superposition holds. 

2. Uniqueness : the solution is uniquely determined for t > 0 by the explicit formula 

/ 0 ° 1 2 

Q(* - y, t)f(y)dy, Q{z ) = -^= e ^r > o. 

v47t at 

3. Existence for large enough set of addmisible initial data: bounded initial data /(x) 
can be prescribed (or at least j/(x)| < e ai2 ), and the corresponding solution is C°° - 
in fact , u(x, t > 0) is analytic because of exponential decay in Fourier space. 

4- The maximum principle : follows directly from the representation of u(x, t) as a convo- 
lution of f(x) with the unit mass positive kernel Q(z). 

5. Energy decay : as usual we may proceed in one of two ways. 

The Fourier Method . We start with 

II Oil, < 2iry; \!(k)f . Max*[|Jf . < ConsU- - 

(since the quantity inside the above brackets is maximized at k 2 at - s ). The last a priori 
estimate shows that the parabolic solution becomes infinetly smoother than its initial data 
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is as we "gain” infintely many s-derivetives, and at the same time, higher derivetives decay 
faster as t f oo. Alternatively, we can work with the 

Energy Method. 


Multiply (2.2.1) by u and integrate to get 


( 2 . 2 . 6 ) 


1 ^ || it 2 

2^NIl 2 (*) 


< — a\\u x \ 


2 

M x ) 


and in general 
(2.2.7) 


1 d d 3 u 

2 dt dx s 


2 

l 2 


< —Const. 


<9 5+1 u 

dx 


2 

l 2 


successive integration of (2.2.7) yields (2.2.5). 

Turning to general case, we consdier mth-order P.D.E.’s of the form, 


( 2 . 2 . 8 ) 


du 

dt 


P(x, t, D)u, P(x, t,D)=J2 0^'- 

lil=o 


We say that the system (2.2.8) is weakly parabolic of order a if 
(2.2.9) ||^ 7 u(x,i)|| L2 < Const.r |i|/ “||u(x,0)|| L2(x) . 

In the constant coefficients case this leads to the Garding-Petrovski characterization of 
parabolicity of order /?, requiring 


ReA 


m 

P(ik) = £ A^ky 


lil=° 


< -A-\kf + B. 


Remark : Generically we have a = /3 = m the order of dissipation which is necessarily 
even. 

The extension to the variable coefficients case (with Lipschitz continuous coeeficients) 
may proceed in one of two ways. Either, we freeze the coeeficients and apply the Fourier 
method to the constant coeficients problems; or, we may use the energy method, e.g., inte- 
gration by parts shows that for 



with Aj + Aj > 8 > 0, and Bj = Bj , the corresponding systems (2.2.8) is parabolic of 
order 2. 

Example: u t = au xx + u xxx is weakly parabolic of order two, yet does not satisfy Petrovski 
parabolicity. 
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2.3. Well-Posed Time-Dependent Problems. 


Hyperbolic and parabolic equations are the two most important examples of time-dependent 
problems whose evolution process is well-posed. Thus, consider the initial value problem 

du 


(2.3.1) 


dt 


= P(x } t , D)u. 


We assume that a large enough class of admissible initial data 
(2.3.2) u(x y t = 0) = f{x) 

there exists a unique solution, u(x,t). This defines a solution operator, E{i y r) which de- 
scribes the evolution of the problem 


(2.3.3) 


u(i) = E{t y t)u(t). 


Hoping to compute such solutions, we need that the solutions will depend continuously in 
their initial data, i.e., 


(2.3.4) 


|u(t) — v(i)|| < Const7||u(0) — v(0)||h* 0 < t < T. 


In view of linearity, this amounts to having the a priori estimate (boundedness) 

(2.3.5) ||ii(t) = E(t y r)u(r)\\ < Constr||u(r)||/f-, 0 < t < T, 
which includes the hyperbolic and parabolic cases. 

Counterexample: (Hadamard) By Cauchy-Kowaleski, the system 

,A = 

has a unique solution for arbitrary analytic data, at least for sufficiently small time. Yet, 
with initial data 

(2.3.6) 

we obtain the solution 


du du 

- + A^ = 0, u = 


Ui 

U 2 


0 +1 

-1 0 


u i{ x > 0) = — , u 2 (x,0) = 0, 


n 


cosh nt sin nx sinhntcosnx 

Ui(x y t) = , u 2 {x i t) = — 


n 


n 


which tends to infinity ||tx( • , <)|| n _ 


oo, while the initial data tend to zero. Thus this 


Laplace equation, -ft- + = 0, is not well-posed as an initial- value problem. 

Finally, we note that a well-posed problem is stable against perturbations of inhomogeneous 
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data in view of the following 

Duhamel’s principle . The solution of the inhomogeneous problem 


(2-3.8) — = P(x,t,D)u + F(x,t) 

is given by 

(2.3.9) u(t) = E(t,0)u(0) + I' E(t,r)F(T)dT. 

J T= 0 

Proof : We have 

~u(t) = j t m,0)u(0)] + j t [J^E(i,T)F(t)dT\ 

= P(x,t, D){E(t,0)u(0)] + E(t,t)F(t) + /' ^\E(i,r)F(t)\dr 

Jt = 0 Ot 

= P{x,t,D)[E{t,0)u(0) + f E(t,T)F(T)dT] + F{t) = P(x,t,D)u{t) + F. 

J T— 0 

This implies the a priori estimate 

(2.3.10) \\u(t)\\ < Const r ||u(0)|| + Const T f* ||F(r)|| w .rfr > 0 < t < T, 

J T— 0 

as asserted. 
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3. THE FOURIER METHOD FOR HYPERBOLIC AND PARABOLIC EQUA- 
TIONS 


3.1. The Spectral Approximation 

We begin with the simplest hyperbolic equation — the scalar constant-coefficients wave 
equation 


(311) 

subject to initial conditions 
(3.1.2) 


du du 
dt dx 


u(x,0) = /(x), 


and periodic boundary conditions. 

This Cauchy problem can be solved by the Fourier method: with /(x ) 5^ — 00 /(*)«“■ 

we obtain after integration of (3.1.1), 


g 

— u(k,t ) = ikau(k,t), 
ot 


u(k,t) = e ikat f(k), 


(3.1.3) 

with solution 

(3.1.4) 
and hence 

(3.1.5) u{x,t) = '£e ikat f(k)e ikx = £ /(fc)e ifc(x+at) = f(x + at). 

k k 

Thus the solution operation in this case amounts to a simple translation 

E(t, t)u(x, t) = u(x + a(t- t ), t), \\E(t, t)|| = 1. 

This is reflected in the Fourier space, see (3.1.4), where each of the Fourier coefficients has 
the same change in phase and no change in amplitude; in particular, therefore, we have the 
a priori energy bound (conservation) 


(3.1.6) 


||«(-,oil 3 = 2*£ M*. <)l 2 = 2 * £ |/(t)| ! = ||/(- 


We want to solve this equation by the spectral Fourier method. To this end we shall ap- 
proximate the spectral Fourier projection of the exact solution Sun = S Nu[x,t). Projecting 
the equation (3.1.1) into the N-space we have 


(3.1.7) 


•N 


du 

dt 


= S 


N 


\a 


du 

dx 
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Since 5# commute with multiplication by a constant and with differentiation we can write 
this as 

/ q 1 q \ Ollff 

(3 ' L8) ~dt = a ~dx' 

Thus un = Snu satisfies the same equation the exact solution does, subject to the approxi- 
mate initial data 


(3-1-9) u N {t = 0) = S N f. 

The resulting equations amount to 2N + 1 ordinary differential equations (O.D.E.) for the 
amplitudes of the projected solution 

(3-1-10) —u N (k, t) = ikau N (k, t ), — N < k < N, 

subject to initial conditions 


(3-1-11) tW(M) = /(*). 

Since these equations are independent of each other, we can solve them directly, obtaining 
(3-1-12) ii N (k, t) = e ikat f(k) 


and our approximate solution takes the form 


(3.1.13) 


N 

Uff(x,t)= ^2 f(k)e lk( - x+at \ 


k=-N 

Hence, the approximate solution u N (x, t ) = f N ( x + at) satisfies 
(3-1-14) u(x, t ) - u N (x,i) = E(t, 0 )/(x) - E(t, 0 )S N f(x) 


and therefore, it converges spectrally to the exact solution, compares (1.1.26), 


(3.1.15) 


K*)-M*)ll < W, m-s N )f(x)\\< 

< ||(r- S N )f(x)W < Constll/lln. • -b.. 


Similar estimates holds for higher Sobolev norms; in fact if the initial data is analytic then 
the convergence rate is exponential. In this case the only source of error comes from the 
initial data, that is we have the error equation 


(3.1.16) 


—[u - Uiv] = - ujv] 
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subject to initial error 
(3.1.17) 


u — — 0 ) — / — /at- 


Consequently, the a priori estimate of this constant coefficient wave equation 


(3.1.18) 


j|u - u N (t)\\ < Constrll/ - fs\\ < Const.il/llH* • ^7 Constr = L 


Now let us turn to the scalar equation with variable coefficients 
(3.1.19) Y t = a ( x , °( x > t) = 2*- periodic. 

This hyperbolic equation is well-posed: by the energy method we have 

=- f u x a(x,t)u-f a z (x,t) u 2 


(3.1.20) 
and hence 
(3.1.22a) 
with 
(3.1.226) 


-— [ u\x,t)dx= f ua(x,t)u x = J a x (x,t)u 2 (x,t)dx, 

2 dt J x ** 


||ti(*,t)||L,(*) < Constr • ||/(x)|| 


Consty = e MT , M = Max x> t[— a x (x, i)]. 


In other words, we have for the solution operator 


\\S(t, t)u(t) || M x) < e M(t t) ||w(t)||l 2 (*) 


and similarly for higher norms. As before, we want to solve this equation by the spectral 
Fourier method. We consider the spectral Fourier projection of the exact solution u N = 
Sfju(x,t)] projecting the equation (3.1.19) we get 

d q 

(3.1.23) ^7 un = 


a(x , t) 


du 

dx 


Unlike the previous constant coefficients case, now Sn does not commute with multipli- 
cation by a(x,t), that is, for arbitrary smooth function p(x,t) we have (suppressing tqne 

dependence) 


(3.1.24a) 


S N a(x)p(x) 


E ( E a(fc-i)Ki))^ 

k=-N \j--00 } 
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while 


oo / N 

(3.1.24 b) a(x)S N p{x)= £ E K k ~j)pU) 

fc— — oo \j= — N 

Thus, if we exchange the order of operations we arrive at 
(3.1.25) ^ = a(x, t)^-~ [a(z, t)S N - S N a(x, t)}^. 

While the second term on the right is not zero, this commutater between multiplication and 
Fourier projection is spectrally small, i.e., 

||Sjva(x)p(x) - a(x)5' J vp(a:)||L 2 (x) = 



(3.1.26) 


\\(Sn ~ 7)a(x)p(x) + a( x ){I ~ < 

< Const. ||a(x)p(x)|| H . • + Const. ||a(x)|| Loo(l) • ||p(z)||tf' • jj~ 


and so we intend to neglect this spectrally small contribution and to set as an approximate 
model equation for the Fourier projection of u (x,t) 

d y N / ,n dvif 

Yet, the second term may lie outside the N-space, and so we need to project it back, thus 
arriving at our final form for the spectral Fourier approximation of (3.1.19) 

dv N \ 


(3.1.27) 


N 


dv 

~dt 


= S N a(x, t) 


dx 


Again, we commit here a spectrally small deviation from the previous model, for 


(3.1.28) 


1 


||(7 - S w )ap(z)||L 2 (x) < Const||a(x)p(x)|| H . ■ jjj 


The Fourier projection of the exact solution does not satisfy (3.1.22), but rather a near by 
equation, 


(3.1.29) 


dux 

~dT 


S N (a(x,t)^]+F N (x,t) 


where the truncation error . F^{x t t ) is given by 
(3.1.30) F N (x,t) = S N 


du 

a(x, t)(I — 
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The truncation error is the amount by which the (projection of) the exact solution misses 
our approximate mode (3.1.27); in this case it is spectrally small by the errors committed in 
(3.1.26) and (3.1.18). More precisely we have 

(3.1.31) \\F N (X } t)\\L 2 (x) < \\ a ( X j 0H^2(*) ' IM|h* + 1 

depending on the degree of smoothness of the exact solution. We note that by hyperbolicity, 
the later is exactly the degree of smoothness of the initial data, i.e., by the hyperbolic 
differential energy estimate 

(3.1.32) ||/ J Vr(*,0IUa(*) - ll a ( x >*)IUa(*) ' ll/lltf*+i ' jj 7 

and in the particular case of analytic initial data, the truncation error is exponentially small. 

From this point of view, our spectral approximation (3.1.27) satisfies an evolution model 
which is spectrally away from that of the Fourier projection of the exact solution (3.1.29). 
This is in addition to the spectrally small error we commit initially, as we had before 

(3.1.33) uj/(t = 0) = Sitf = 

We now raise the question of convergence. That is, whether the accumulation of spectrally 
small errors while integrating (3.1.27) rather than (3.1.29), give rise to an approximate 
solution ujv(i, t ) which is only spectrally away from the exact projection ujv(z, t). We already 
know that the distance between u^(x, t ) and the exact solution u(x, t) — due to the spectrally 
small initial error - is spectrally small as we have seen in the previous constant coefficient 
case. 

To answer this convergence question we have to require the stability of the approximate 
model (3.1.27). That is, we say that the approximation (3.1.27) is stable if it satisfies an a 
priori energy estimate analogous to the one we have for the differential equation 

(3.1.34) || y w(*)|| < Const.e Mt ||vjv(0)||. 

Clearly, such a stability estimate is necessary in any computational model. Otherwise, the 
evolution model does not depend continuously on the (initial) data, and small rounding 
errors can render the computed solution useless. And on the positive side we will show that 
the stability implies the spectral convergence of an approximate solution Indeed 

the error equation for = u^(t) — vpf(t) takes the form 

(3.1.35) *^L = S n a(x,t)^¥- + F N (x,t) . 

3 We note that in the previous constant coefficient case, the approximate model coincides with the differ- 
ential case, hence the stability estimate was nothing but the a priori estimate for the differential equation 
itself. 
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Let Eff(t,r) denote the evolution operator solution associated with our approximate model. 
By the stability estimate (3.1.34) 

(3.1.36) ||Ej V (t,r)u^(r)|| < Conste M(t-r) ||u JV (r)||. 


Hence, by (3.1.36) together with DuhammePs principle we get for the inhomogeneous error 
equation (3.1.35) 


(3.1.37a) 

and 


e^(t) = E^(t ) 0)e^(0) + f t)Fn(t)<1i 

J r= 0 


(3.1.376) ||ejv(t)|| < Const. e Mt ||ejv(0)||L 2 (x) + / H-FaK®, t )||lj(*)^ • 

J T— 0 

In our case e^O) — f N — Sfw = 0, and the truncation error Fn(x,t) is spectrally small; 
hence 


(3.1.38) IJe/v = u/y(i) — ujv(i)|| < Const. e Mt • 

where the constant depends on ||a(x, f)||^ QO (!) and ||/|| i.e., restricted solely by the 
smoothness of the data. In the particular case of analytic data we have exponential conver- 
gence 


(3.1.39) ||ejv'(f) = u N (t) — u^(f)|| < Const. e Mt • e . 


Adding to this the error between Utf(t) and u(t) (- which is due to the spectrally small error 
in the initial data between f n and /) we end up with 

jjr for H J+1 initial data 

e~ nN for analytic initial data 

To summarize, we have shown that our spectral Fourier approximation converges spectrally 
to the exact solution, provided the approximation (3.1.27) is stable. 

Is the approximation (3.1.27) stable? That is, do we have the a priori estimate (3.1.34)? 
To show this we try to follow the steps that lead to the analogue estimate in the differential 
case, compare (3.1.20). Thus, we multiply (3.1.27) by v^(x } t) and integrate over the in- 
period, obtaining 

(3.1.41) ~ f"V 2 N (x,t)dx = + J*v N (x,t)S N ^a(x,t)^Ej dx. 


(3.1.40) ||u(t) - ujv(i)|| < Const. e Mt • { 
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But is orthogonal to (7 — S n) |<x(x, t)^j- L J so adding this to the right-hand side of 

(3.1.41) we arrive at 

( 3 .L 42 ) “ J* v 2 N (x,t) = J* v N (x, t)a(x, i)^-dx 

and we continue precisely as before to conclude, similarly to (3.1.22), that the stability 
estimate (3.1.34) holds 

(3.1.43) ||vjv(t)|| ^ Const. e M *||u^(0)||, M = Max X|t [— a x (sc, t)}. 


In the constant coefficient case the Fourier method amount to a system of (2 N + 1) 
uncoupled O.D.E.’s for the Fourier coefficients of v jy = u jy which were integrated explicitly. 
Let’s see what is the case with problems having variable coefficients say, for simplicity, 
a = a(x). Fourier transform (3.1.22) we obtain for t — f)jv(fc, <) - the kth-Fourier 
coefficient of v^(x y t) = J2k=-N t)e tkx , 


(3.1.44) dv{k f t) _ a(Ai - j)ijv(j, t), -N < k < N. 

dt j=-N 

In this case we have a (2 N + 1) X (2 N + 1) coupled system of O.D.E.’s written in the 
matrix- vector form, consult (1.2.46) 


(3.1.45) — u(t) = .AAu(t), u(t) = 

dt 


v(N,t) J 


A kj = d(k - j), A = diag(ifc). 


We can solve this system explicitly (since a (•) was assumed not to depend on time) 
(3.1.46) v(t) = e AAt t>(0); 


that is, we obtain an explicit representation of the solution operator 
(3.1.47) E N (t,r) = 7^V a (‘- t )7V, A - A N ,A = A N 


where Fn denote the spectral Fourier projection 


(3.1.48) 


FnVn(x) 


v(-N) 

v(N) 


We note that in view of Parseval identity ||i'Vu,iv (®)||2 = ||ujv(x)||£ J ( x ) (modulo factorization 
factor), hence, stability amounts to having the a priori estimate on the discrete symbol 
EN(t, T ) — e^ NK ^~ T \ requiring 


(3.1.49) 


e i„A(t-r)|| < Const.e M(t_T) . 
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The essential point of stability here, lies in having a uniform bound in the RHS of (3.1.49) 
independent on the order of the system, e.g., a straightforward estimeate of the form 


(3.1.50) 


r)|| <- gll^jvIHIAIKt— r) 


will not do because || A/sr || oo. The essence of the a priori estimate we obtained in 
(3.1.22), and likewise in (3.1.42), was that the (unbounded) operator P(x,t,D) = a{x,t)d x 
is semi-bounded , i.e., 


Re 


a(x, t ) 


dx 


. . d d . . . . 

a ^ X,t ^dx ~ dx^ X,t ^ 


1 / X 

— 2 <Ix \ a '> )i 


namely, compare (2.1.26) 


r r \ 

Re a(x,t)— u,uj < M||u||^ (x) 
^ L J / l 2 (x) 


and likewise for Re j. In the present form this is expressed by the sharper 

estimate of the matrix exponent, 4 compare (3.1.50) 


(3.1.51) 


| e ii N A(t-r)|| < g| j | -(t—T) . 


This time, (| Re^4 || like the R e[P(x,t> D)], is bounded. Indeed, [ReAAjfcj = |[a(fc — j)ij + 
a[j — k)ik], and since a[x } t) is real (hyperbolicity!) then a(p) = a(— p), i.e., 


(3.1.52) 


[ReAA] kj = -i(j - k)a(k - j) - N < j, k < N. 


Thus, ReAA is a Toeplitz matrix, namely its (A:, j) entry depends solely on its distance from 
the main diagonal k — j\ we leave it as an exercise - using our previous study on circulent 
matrices in (1.2.43) - to see that its norm does not exceed the sum of absolute values along 
the, say, zeroth (j = 0) row, i.e., 

1 N 

(3.1.53) n ~ ' J 


HReA»A|| < - £ \ka(k)\ 


k=N 


which is bounded, uniformly with respect to A/ , provided a(x } t) is sufficiently smooth, e.g,, 
we can take the exponent M to be 


M 


(3.1.54) 


= 5 £ M*)l < IJtJ' KW- E 1< 

" ^ * r xr A 


N 

c 

k=-N 


7 r 


— ' || a xx( X ) 0 ||l 2 (x) 


4 To see this, use Duhammers principle for ^ — ReAAu(i) + F(t) where F(t) = HmAAe AAt or integrate 
directly. 
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which is only slightly worse than what we obtained in (3.1.43). 

A similar analysis shows the convergence of the spectral-Fourier method for hyperbolic 
systems. For example, consider the N x N symmetric hyperbolic problem 


(3.1.55a) 


du .du . . . 

— = A(x,t)- — \- B{x,t)u, with symmetric A{x,t). 

C/£ \j I 


We note that if the system is not in this symmetric form, then (in the 1-D case) we can bring 
it to the symmetric form by a change of variables, i.e., the existence of a smooth symmetric 
H(x, t) such that H(x, t)A(x, t ) is symmetric, implies that for w(x, t ) = T -1 (x, t)u(x, t ) with 
H — (J'- 1 )*J 1 - 1 we have, compare (2.1.15b) 

3%d 3xu 

(3.1.556) = T -1 (x, t)A(x, t)T(x, t + C(x, t)w(x, t) 

where T _1 (x, t)A(x, t)T(x, t) = T*(x,t)H(x,t)A(x,t)T(x,t) is symmetric, and C(x,t ) = 
B(x,t)+ aT dt --(x, t)—T~ l (x, t)A(x , t)|^(x, f). The spectral Fourier approximation of (3.1.55a) 
takes the form 


(3.1.56) 


dv 


'N 


dt 


= S N A(x, t) 


du 


N 


dx 


+ SNB(x,t)v N (x,t). 


Its stability follows from integration by parts, for by orthogonality 


(3.1.57a) — — J v 2 N (x,t)dx = J t >^A(x, t) -^- dx + J u^B(x,t)uffdx <mJ v^(x,t)da 

dA(x, t ) 


where 
(3.1.576) 
and hence 
(3.1.58) 


M = Max xt 


dx 


+ Re5(x, t ) 


IMOII^W < e"*llw(0)||. 


The approximation (3.1.56) is spectrally accurate with (3.1.55) and hence spectral conver- 
gence follows. The solution of (3.1.56) is carried out in the Fourier space, and takes the 
form 

d N 

(3.1.59) 3tv(M) = M k ~ ~N <k< N, 

at j=-N 

which form a coupled (2 N + 1) X (2 N + 1) system of O.D.E.’s for the (2iV + l)-vectors of 
Fourier coefficients v(k,t). 

There are two difficulties in carrying out the calculation with the spectral Fourier method. 
First, is the time integration of (3.1.59); even in the constant coefficient case, it requires to 
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compute the exponent e AAt which is expensive, and in the time-dependent case we must 
appeal to approximate numerical methods for time integration. Second, to compute the 
RHS of (3.1.59) we need to multiply an (2 N + 1) x (2 N + 1) matrix, AA by the Fourier 
coefficient vector which requires 0(N 2 ) operations. Indeed, since A is a Toeplitz matrix and 
A is diagonal, we can still carry out this multiplication efficiently, i.e., using two FFT’s which 
requires O(iVlogiV) operations. Yet, it still necessitates to carry out the calculation in the 
Fourier space. We can overcome the last difficulty with the pseudospectral Fourier method. 

Before leaving the spectral method, we note that its spectral convergence equally applies 
to any P.D.E. 

du 

(3.1.60) = P(x,t,D)u 

with semi-bounded operator P(x y i, D), e.g., the symmetric hyperbolic as well as the parabolic 
operators. Indeed, the spectral approximation of (3.1.60) reads 

dv 

(3.1.61) -7^ = S N P( x, t, D)v n . 

Multiply by v n and integrate - by orthogonality and semi-boundedness we have 

(3.1.62) f v 2 N (x y t)dx = Re(u^, P(x } t, D)v N ) < M [ v 2 N (x>t)dx. 

2 dt Jx Jm 

Hence stability follows and the method converges spectrally. 
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3.2. The Pseudospectral Approximation 


We return to the scalar constant coefficient case 

du du 


dt a dx 


(3.2.1) 

subject to periodic boundary conditions and prescribed initial data 

(3.2.2) w(x,0) = /(x). 

To solve this problem by the pseudospectral Fourier method, we proceed as before, this time 
projecting (3.2.1) with the pseudospectral projection tpN> to obtain for u n = i p^u(x,t) 


(3.2.3) 


dux . ( du\ 

= Vn a— . 


dt 


dx 


Here, ipN commutes with multiplication by a constant, but unlike the spectral case, it does 
not commute with differentiation, i.e., by the aliasing relation (1.2.2) we have 


(3.2.4a) 
where as 
(3.2.46) 


£ (kp{k))e ikx = £ 52i[k + j(2N+l)]p[k + j(2N + l)]e 

° X kzz~N k=~N j 


ikx 


N 


N 


£-+»P- E **E +!)]<"■. 

OX k— — N k= — N j 

The difference between these two operators is a pure aliasing error, i.e., we have for = 
Sn + An, see (1.2.13) 


*4 - 4 (vw) 


' d 
An, — 

dx 


N 


P= £ E*[ fc +i ( 2 N + l)]p(k + j(2N+l)]e 

k=-N j^O 


ikx 


which is spectrally small. Sacrificing such spectrally small errors, we are led to the pseu- 
dospectral approximation of (3.2.1) 


(3.2.5) 

subject to initial conditions 

(3.2.6) 


dvN dv?{ 
= a- 


dt 


dt 


v N (t = 0) = ip N f. 


Here, Vn — VN(x,t ) is an N-degree trigonometric polynomial which satisfies a nearby equa- 
tion satisfied by the interpolant of the exact solution TpNu(x,t). That is, un = ipNu{x,t) 
satisfies (3.2.5) modulo spectrally small truncation error 


(3.2.7) 


~W ~ <1 ~^x' + FN ( X,t }’ F N(x,t) = aipN 
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where by (3.2.3), Fff(x,t) 
small 


= a 


j, and by (1.2.17) it is indeed spectrally 


( 3 - 2 - 8 ) ll-PW(*,OII < l a l II^P- V’wHIl] < M |M| h .+i-T 

The stability proof of (3.2.5) follows along the lines of the spectral stability, and spectral 
convergence follows using Duhammel’s principle for the stable numerical solution operator. 
That is, the error equation for e n = is 

whose solution is 


(3.2.10) e N (t) = E N (t,0)(f N - ip N f) + f E N {t,T)F N {x,r)dT. 

J T= 0 

Hence, by stability 

(3.2.11) || e w(t)|| < Const. • \\u\\ H .» ~ < Const.e m \\f\\ H ^ • T- ; 
this together with the estimate of the pseudospectral projection yields 

{ ^7 for H 3+1 initial data 
e~ nN for analytic initial data 

To carry out the calculation of (3.2.5) we can compute the discrete Fourier coefficients v(k, t ) 
which obey the O.D.E., 


(3.2.13) —(k,v) = ikav(k,t), 

at 

as was done with the spectral case; alternatively, we can realize our approximate interpolant 
djv(x, t ) at the 2iV -f 1 equidistant points x v = vh, and (3.2.5) amounts to a coupled ( 2N + 1) 
- O.D.E. system in the real space 


(3.2.14) 


dv N dv N 

-{x V) t) = a-^— (x 


dt 


x„, t ) v = 0, 1, • • • ,2 N. 


(3.2.15) 


u w (x^,0) = f(x„). 


Let us turn to the variable coefficient case, 

du 


(3.2.16) 


/ \ 

at =o(l '%' 
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The pseudospectral approximation takes the form 

dv N 


(3.2.17a) 


dt 


= i’N 


a(x, t ) 


9vn 

dx 


(3.2.17 b) 


V N (x v , 0) = /(£„). 


It can be solved as a coupled O.D.E. system in the Fourier space, but just as easily can be 
realized at the 2 N -f- 1 so-called collocation points 


(3.2.18a) 


dv N {x v , t) 


. dvt 


dt Q-{x vy f) (® — x V) t) 


(3.2.186) 


v N {x v , t = 0) = /(x„). 


The truncation error of this model is spectrally small in the sense that u n = satisfies 




dt 


(3.2.19) 
where 

(3.2.20) Fn(x, t) = ipx 

is spectrally small 


4>N 


a(x, t ) 


du 


N 


dx 


+ F N (x,t ) 


. . du 

a(l '% 


V’jv 


Q 

a (x,t)—(ip N u) 


II^OM)!! ^ II i>N 


(3.2.21) 


~ M u \\\ < 


< || a (^,i)||l,“ ' ||/||tf*+i ' 


Hence, if the approximation (3.2.12) stable, spectral convergence follows. Is the approxima- 
tion (3.2.12) stable? The presence of aliasing errors is responsible to a considerable difficulty 
in proving this, and currently this is an open question. If we try to follow the differential 
and spectral recipe, we should multiply by Vflf(x,t) and integrate by parts. However, here 
U//(cc,t) is not orthogonal to (7 — V’w)[‘ ' '] which otherwise would enable us to follow the 
differential estimate of / v^(x,t)a[x,t)^^-(x,t)dx in terms of f x v%(x ,t)dx\ more precisely, 
we have for ip N = S N -f A N that I - ip N = I - S N - A N where / v N (I - Sn)[- ■ -}dx = 0, 
compare (3.1.41), (3.1.42); yet / -jdi leaves us with an additional contribution which 

cannot be bounded in terms of J x v 2 n (x , t)dx in order to end up with the stability proof with 
Gronwall’s inequality. To shed a different light on this difficulty, we can turn to the Fourier 
space; we write (3.2.17) in the form 


(3.2.22) 


, .dv N 

IT = a(l) -57 
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and Fourier transform to get for the kth Fourier coefficient 


(3.2.23a) 



N 

E a(k - j,t)ijv(j,t) 


j=-N 


i.e., 

(3.2.236) 


^-v(t) = A N Av(t) A kj = Y & i k -3 + P( 2N + !)]• 
at p 


This time, ReA^A is unbounded. This difficulty appears when we confine ourselves to the 
discrete framework: multiplying (3.2.18a) by u(x v ,t) and trying to sum by parts we arrive 
at 


(3.2.24) 


1 d dv 

— — ^ ^(x^, t) = ^ a{x U} t)v(x V} t ) 0 


E — 

^dx 


-a(x,*)u 2 (x,i) 


n 0, { X V> 0 U Af( a '»'> 

1 J ^ 


but the first term on the right does not vanish in this case - it equals, by the aliasing relation, 
to 


E — 

V dx 


|i a(x,t)v 2 (x,t ) 


= / ^[• • ] + Ep ( 2JV + 1 )5«>[P'(2W + 1)] 

x~ x i/ P5^0 


and a loss of one derivative is reflected by the factor 2iV + l inside the right summation. There 
are two main approaches to enforce stability at this point: skew-symmetric differencing and 
smoothing. We discuss these approaches in the next two subsections. 
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3.3 Skew-Symmetric Differencing The essential argument of well-posedness for symmetric 


hyperbolic systems with constant coefficients is the fact that (say in the 1-D case) P(D) = 

is a skew-adjoint operator what is loosely called “integration by parts” With variable 

coefficients this is also true, modulo low-order bounded terms, i.e., 

d 1 " Q Q | 

(3.3.1) P(x , t, D ) = A(x, "k ~Qx^ a ^ X} ^ ^ ~~ 

The stability proofs of spectral methods follow the same line, i.e., we have in the Fourier 
space, compare (3.1.45), (3.1.52) 

(3.3.2) A n A = - [A/vA — A Avr] + ^ [AatA + AMjv] 

and stability amounts to show that the second term in (3.3.2) is bounded: for then we have 
in (3.3.2) (as in (3.3.1)) a skew-adjoint term with an additional bounded operator. The 
difficulty with the stability of pseudo-spectral methods arises from the fact that the second 
term on the right of (3.3.2) is unbounded, 

(3.3.3) lim^ooll^Az/A + A* A N )\\ T °°- 

To overcome this difficulty, we can discretize the symmetric hyperbolic system (again, say 
the 1-D case) 

/ \ s du 

(3.3.4a) — = A(x, i ) — 

when the spatial operator is already put in the “right” skew-adjoint form, compare (3.3.1), 

(3.3.46) 1ft A ( X,t ^ + J^UOMh) ~^A x (x,t)u. 

The pseudospectral approximation takes the form 

(3.3.5) = 1 il> N + ^-Tf> N {A{x,t)u N ) - )^ n {A x {x, t)u N ). 

In the Fourier space, this gives us 

(3.3.6) §=I[i„A + Ai N ]v-i^%. 

Now, Ax A + AAjf is symmetric because A is, is bounded and stability follows. 
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3.4. Smoothing 


We have already met the process of smoothing in connection with the heat equation: 
starting with bounded initial data, /(x), the solution of the heat equation (2.2.1) 

1 

(3.4.1) u(x, t) = Q * /(x), Q(x) = —== e , t> 0 

\/47 ra 

represents the effect of smoothing /(x), so that u(*,f > 0)eC°° (in fact analytic) and u(x,i j 
0) = f(x). 

A general process of smoothing can be accomplished by convolution with appropriate 
smoothing kernel Q(x) 

(3-4.2) f e (x) = Q.(x) * f{x) 

such that 

(3.4.3a) Q c {x) * f{x) 


is sufficiently smoother than /(x) is, and 

(3.4.36) Q t {x)*f(x )— * f(x). 


tf— ►0 


With the heat kernel, the role of e was played by time t > 0. A standard way to construct 
such smoothers is the following. We start with a C’-function supported on, say, (-1,1), such 
that it has a unit mass 


J Q(x)dx = 


1 


(3.4.4a) 

and zero first r moments 

(3.4.46) J x*<j>(x)dx = 0, j = 

Then we set Q c (x) = \Q (j') and consider 


, r. 


(3.4.5) 


fe(x) = Q c (x) * f(x), e > 0. 


Now, assume / is (r + 1) - differentiable in the e neighborhood of x\ then, since Q c (x ) is 
supported on (—e,s) and satisfies (3.3.4a) as well, we have by Taylor expansion 


f(x) - Q c {x ) * f(x) = f Q c (y)[f(x) - f(x - y)]dy = 

J \y\<* 


(3.4.6) 


= L Q - (y) \i t f fW(x)+ v- r9 r+,(0 


dy. 
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The first r moments of Q c {y ) vanish and we are left with 

(3.4.7) |/(x) - Q e (x) * /(x) | < Const. Max| y _ x |< e |/ (r+1) (i/) • £ r+1 , 

i.e., f e (x) converges to f(x ) with order r + 1 as e — * 0. Morevoer, f e ( x) is as smooth as (f)(x ) 
is, since 

(3.4.8) hi*) = J y -Q (^ 7 ^) f^ dy 

has many bounded derivatives as Q has, i.e., starting with differentiable function / of order 
7 * -)- 1 in the neighborhood of x, we end up with regularized function / e (x) in C , s > t. 


Example: For C°° regularization - choose a unit mass C°° kernel, see Figure 6, 

Qoe~^-* 7 , Jx | < 1 , 

with Q 0 such that / <^>(x)dx = 1. 

0, |x| > 1 

Then /«(x) = Q e (x ) * /(x) is a C°° regularization of /(x) with first order convergence rate 


(3.4.9a) 


4 >{ x ) = 



Figure 6: 


(3.4.9&) |/(x) - / £ (x)| < Const. Maxi^^^el/^y)! • e -* 0. 

To increase the order of convergence, we require more vanishing moments which yield more 
oscillatory kernels as in Figure 7. 
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Figure 7: 

We note that this smoothing process is purely local - it involves e-neighboring values of 
C r+1 function /, in order to yield a (^'-regularized function f e (x ) with f e (x ) — >f(x). The 
convergence rate here is r + 1 . 

We can also achieve local regularization with spectral convergence. To this end set 
(3.4.10a) Qn{%) = p(x)D N { x) 

where p(x) is a (^“-function supported on ( — £,£) such that, see Figure 8, 

(3.4.106) p ( 0) = 1. 



Figure 8: 
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Consider now: 


(3.4.11) 


f N (x) = Qn * /(*)■ 


Then /at(x) is C°° because Q N is; and the convergence rate is spectral, since by (1.1.25) 


(3.4.12a) 


f(x) -Q n * f(x) = f(x) - f D N (y) P (y)f{x - y)dy 

c/|y|<c 

= f{x) - p(y)f(x - y)| v =o + residual 


where 

(3.4.126) | residual | < Const. \\p(y)f{x - y)|| *•(-*.«) for an ^ 3 > °- 

Thus, the convergence rate is as fast as the local smoothness of / permits; of course with 
p= 1 we obtain the C^-regularization due to the spectral projection. The role of p was to 
localize this process of spectral smoothing. 

We can as easily implement such smoothing in the Fourier space: For example, with the 
heat kernel we have 


(3.4.13) u(k,t) = e~ akH f(k ) 

so that u(k,t) for any t > 0 decay faster than exponential and hence u(x,i > 0) belong 
to H * for any s (and by Sobolev imbeddings, therefore, is in C°° and in fact analytic). In 
general we apply, 

°o 

(3.4.14) fc(x) = X] Q e {k)f(k)e' kx 

k— — oo 


such that for f e (x) to be in H 3 we require 

°o A 

(3.4.15a) (1 + |fc| 2 m<WI 2 (/(*)l — Const. 

kxx — a. 

and r + 1 order of convergence follows with 

(3.4.156) 1 4>c(k) - 1| < Const. (eA:) r+1 . 


Indeed, (3.4.15b) implies 


(3.4.16) 


\f(x)-fe(x)\ < 

< 


OO 

Const. £ r+1 £ \k T+1 f(k)e' kx \ 


k= — at 

Const. Max|/ (r+1) | • e r+1 . 
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Note: Since 4> e {k) j 0 we can deal with any unbounded / by splitting £|fc|<|k 0 | + 2|k|>|Jt 0 |- To 
obtain spectral accuracy we may use 


(3.4.17a) 


Q N (k) = 


' si, 

< 

~ smoothly decay to zero 


i*=i < f 

f < |*| < N 


Clearly Qn * f{x ) is C°° and as in Section 1 

(3.4.176) |/(x) - Qn * f(x ) | < \K k ) eikx \ < Const. ||/ 1| • ttttt- 

\k\>N. 

We emphasize that this kind of smoothing in the Fourier space need not be local; rather 
Q c (x) or are negligibly small away from a small interval centered around the origin 

depending on e or jf. (This is due to the uncertainty principle.) 

The smoothed version of the pseudospectral approximation of (3.2.16) reads 

(3.4.18) ^- = -tp N (a(x, t)^{Q *v N )) 

i.e. , in each step we smooth the solution either in the real space (convolution) or in the Fourier 
space (cutting high modes). 5 We claim that this smoothed version is stable hence convergent 
under very mild assumptions on the smoothing kernel Qn( x )- Specifically, (3.4.18) amounts 
in the Fourier space, compare (3.2.3) 

dv 

(3.4.19) ^ = A n AQ n v. 

The real part of the matrix in question is given by 


(3.4.20a) [KcA N AQ N ] kj = i{X k - A,) £ a[fc - j + p(2JV + 1)], - N<k,j<N 

P 

where A Q N = diag A .(zA fc ) 

(3.4.206) i\ k = ikQ N (k) 

is interpreted as the smoothed differentiation operator. Now, looking at (3.4.20a) we note: 


1. For p — 0 we are back at the spectral analysis, compare (3.1.12), (3.1.13) and the real 
part of the matrix in (3.4.20a) — the aliasing free one — is bounded. 

2. We are left with |p| = 1: in the unsmoothed version, these terms were unbounded 
since \X k — Xj\ f co as k j -N or j t N. With the smoothed version, these terms are 
bounded (and stability follows), provided we have 

5 Either one can be carried out efficiently by the FFT. 
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(3.4.21) 


ikQ^(k ) | — * 0 . 

|fc|TW 


For example, consider the smoothing kernel Qn { x ) where 

sin kh . 2 tt 


(3.4.22a) 


Q N (k) — 


kh 


h — 


2N + 1 


This yields the smoothed differentiation symbols 


(3.4.226) 


, . . kh 

At = i sin — 


which corresponds to the second order center differencing in (1.2.36); stability is immediate 
by (3.4.21) for 


_ * 2rrk 

|Afc = 2Ar+1 


2 irk 
2AT+1 


|fc|TJV 


Yet, this kind of smoothing reduces the overall spectral accuracy to a second one; a fourth 
order smoothing will be 


(3.4.23a) 


or 


(3.4.236) 


At = ir 


At = 


. kh . 2kh 

4 sin — sin — — 

h 2 h 


<W*> = s 


6z 


sin 


kh 


Qtf(k) — 


' k ~ h 4 + 2cosfch’ ik\ 

In general, the accuracy is determined by the low modes while stability has to do with 
high ones. To entertain spectral accuracy we may consider smoothing kernels other than 
trigonometric polynomials (= finite difference), but rather, compare (3.4.17) 

fsl, |*| <f 


(3.4.24) 


Q N (k) = 


w smoothly decay to zero y < |A;| < N. 

An increasing portion of the spectrum is differentiated exactly which yields spectral accuracy; 
the highest modes are not amplified because of the smoothing effect in this part of the 
spectrum. 

We close this section noting that if some dissipation is present in the differential model 
to begin with, e.g., with the parabolic equation 


(3.4.25) 


du 

Hi 


d_ 

dx 


a(x, t) 


du 

dx 


a(x, t) > a > 0, 


then stability follows with no extra smoothing. The parabolic dissipation compensates for 
the loss of “one derivative” if first order terms are present. To see this we proceed as follows: 
multiply 


(3.4.26) 


9vn . . d 

-W M = a^ 
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by v^{x u t) and sum to obtain 


(3.4.27) 


2 5 £ = Y^ v N( x u,t)^(a(x l/ ,t)^-(x U) t)). 


Suppressing excessive indices, v N (x v ,t) = v v (t), we have for the RHS of (3.4.27) 


(3.4.28) 



Now, the first sum on the right gives us the usual loss of one derivative and the second 
are compensates with gain of such quantity. Petrowski type stability (gain of derivatives) 
follows. We shall only sketch the details here. Starting with the first term on the right of 

(3.4.28) we have 


(3A29) \ £ k [ a “ v "i£) = \J [aliasing errors] 

while for the second term 


(3.4.30) 


XX(0 


dx 


< 


-a] 


9vn_ 

dx 


1 2 


(*,0 


dx 


and this last term dominates the RHS of (3.4.29). 
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APPENDIX 


A.l. Fourier Collocation with Even Number of Gridpoints 
We assume w(x ) is known at 

(A. 1.1) w„ = w(x„) x v = r + vh v = 0, 1, • • • , 2N 

with h = and 0 < r < h. We use the trapezoidal nodes 


(A1.2) 


1 


w\ 


(*0 = otE'V 

27r u=0 


"... - 


2N 


E 

u=0 


w„e 


— \kx u 


to obtain the pseudospectral approximation 

N 

(A.l. 3) V’ATtw = Y "w(k)e lkx . 

k=-N 


Note : We now have only 2 N pieces of discrete data at the different 2 N grid points 

*o> ^lj ' ' ‘ i I 2 JV-i and they correspond to 2 N waves, as we have a “silent” last mode, i.e., 
with r = 0, k = N, Im[e' kx ] x=Xl/ = isin un = 0. This is a projection, since in view of (A. 1.3) 
tptfW is the interpolant of w(x) at x = x v : 


i> N w(x)\ a 


2 N 


(A.1A) 


N 

= E " 

k=-N L 

2N-1 

= Y w M 


2AT+1 


E ™(x„)e 


— ikx v 


e ikxv 


_ 1 _ 
2 N 


N 


tf ik(p-v)h _ 


U) 


(xj. 


v~0 ’ /c— — N 

The aliasing relation in this case reads - compare (1.2.7) 

oo 

(A. 1.5) U ,(k) = Y e ip2Nr w(k + 2 pN) 

p=-ot 

and spectral convergence follow - compare with (1.2.16) 

(A. 1.6) ||AtfU;(x)||if. < Const, • \\T N w(x)\\ H >, s > 

In the usual sin-cos formulation it takes the form 

(A1.7) 


N 


iprfW — Y, k cos kx + bk sin kx, 


k=0 


<*k 

k 


1 2AA+1 

= 77 E W M 

iV i /=0 


cos kx u 
sin kx u 


0 < k < N. 


Noting that = 0 we have 2 N free parameters ao, {a*, an d to match our data 

at 
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